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Improved Gaussian beam-scattering algorithm ^ 


James A. Lock 


The localized model of the beam-shape coefficients for Gaussian beam- scattering theory by a spherical 
particle provides a great simplification in the numerical implementation of the theory. We derive an 
alternative form for the localized coefficients that is more convenient for computer computations and that 
provides physical insight into the details of the scattering process. We construct a fortran program for 
Gaussian beam scattering with the localized model and compare its computer run time on a personal 
computer with that of a traditional Mie scattering program and with three other published methods for 
computing Gaussian beam scattering. We show that the analytical form of the beam-shape coefficients 
makes evident the fact that the excitation rate of morphology-dependent resonances is greatly enhanced 
for far off-axis incidence of the Gaussian beam. 


1. Introduction 

Although the history of numerical Mie theory compu- 
tations dates back almost to the time of Mie and 
Debye, 1 it was not until a widely published numeri- 
cally stable computer code 2 * 4 could be run quickly on 
low-cost personal computers 5 that plane-wave Mie 
theory computations became commonplace. In re- 
cent years progress has also been made on the 
numerical implementation of other scattering prob- 
lems. For example, a number of computational meth- 
ods for calculation of scattering of a focused Gaussian 
beam by a spherical particle have been devised. In 
these methods the Gaussian beam is expanded either 
in an infinite series of particle waves 6- * or in an 
angular spectrum of plane waves. 9 ” 11 These expan- 
sions have considered not only the Gaussian shape of 
the dominant component of the beam's electric and 
magnetic fields, but also included in an approximate 
way smaller components of the fields induced by the 
variation of the dominant component in the trans- 
verse direction. 1213 

Currently no consensus has been reached as to 
which computational method for Gaussian beam 
scattering is superior to the others or whether one 
method possesses a richness of physical interpreta- 
tion that is not manifest in the others. This is 
because the successes and limitations of each method 
have not yet been fully explored. To this end, in two 
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recent papers 14 - 15 we examined the applicability of 
Gouesbet's localized model 1617 of Gaussian beam 
scattering to the case of tight beam localization. We 
found that for an on-axis beam, i.e., a beam propagat- 
ing along the z axis, the localized model accurately 
describes a focused Gaussian beam with the focal 
waist half-width w 0 satisfying \/{2ttwq) < 0.15. For 
an off-axis beam, i.e., a beam propagating parallel to 
but not along the z axis, the localized model accu- 
rately describes a Gaussian beam for \/(2 t7u;q) < 
0.10, although the accuracy depends to some extent 
on the aspect of the scattering that is being examined. 
In this paper we consider a different application of the 
localized model, namely, the construction of a stable 
and relatively fast-running computer code for Gauss- 
ian beam scattering that can be implemented on a 
personal computer. 

The body of this paper proceeds as follows. In 
Section 2 we give the basic formulas for far-field 
scattering of a focused off-axis Gaussian beam by a 
spherical particle. We also give the localized model 
expressions for the beam-shape coefficients that de- 
scribe the Gaussian beam. We then simplify the 
localized expressions, writing them in terms of either 
Bessel functions or modified Bessel functions of a 
complex argument. In Section 3 we describe algo- 
rithms for numerical computation of the Bessel func- 
tions and other expressions that appear in the far- 
field scattering formulas. In Section 4 we examine 
the computer run time of our Gaussian beam- 
scattering program and compare it with the run time 
of a standard plane- wave Mie theory program. We 
also compare our program with three other computa- 
tional methods for Gaussian beam scattering. Last, 
in Section 5 we show that the analytical expressions 
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for the beam-shape coefficients easily and correctly 
predict the large enhancement in the excitation rate 
of morphology -dependent resonances (MDR’s) by an 
off-axis Gaussian beam focused somewhat beyond the 
edge of a dielectric spherical particle. 


x 





Fig. 1. Focused Gaussian beam that is incident upon a spherical 
particle. The center of the particle is at the origin of the coordi- 
nates. and the center of the beam’s focal waist is at (a) x y, yf \ zf, 
(blxr * 0,yr= 0; and (c)x^ — 0 ,y f * 0. 


2. Localized Model of Gaussian Beam Scattering 

Consider a Gaussian beam focused to the half-width 
w Q at the point (x f ,yfZf) that is incident upon a 
spherical particle whose center is at the origin of 
coordinates. This is shown in Fig. 1(a). The radial 
components of the beam’s electric and magnetic fields 
are E in c rad (r, 0, 4>) and Bin/* 1 (r, 0, <f>), respectively. 
The spherical particle has radius a and refractive 
index n. In the notation of Ref. 18, the far-field 
scattered intensity is 

/(r, 0, <t>) * 2\L<fik 2 r* + 1^®* 4>)| 2 ], ^ 

where the scattering amplitudes Si and S 2 are given 
by 

Si(0, 4>) = 2 2^ 21(1 + 1) 

+ P/m'/ |m| (9)]exp(im(l)) > 

* 1 21 + 1 
S 2 (0, <t>) = 2 2uITT ) 

+ a /m T/ lnit (e)]exp(im<t)). (2) 

In Eqs. (1) and (2) 



is the wave number of the incident beam. The 
partial-wave scattering amplitudes a/ m and p/ m are 
given by 


a im - A lm ai, 

Pim = Bbfih M 

respectively, where CL[ and 6/ are the partial -wave 
scattering amplitudes of plane- wave Mie theory. 
Efficient algorithms for the computation of a { and b\ 
are given in Refs. 3 and 4. The beam-shape coeffi- 
cients A im and B tm are 
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respectively, where ji(kr) are spherical Bessel func- 
tions. The angular functions ir, im| (0) and t^ ot| ( 0) are 
given by 


17,i m| (0) = P/ ml (cos 0), 

* sin 0 

T f |m, (0) = “jP/ !m, (cos 0), (6) 

respectively, where P/ lm| (cos 0) sure associated Legen- 
dre polynomials. 

Calculation of the beam-shape coefficients of Eqs. 
(5) by the use of numerical integration is computation- 
ally the most time-consuming task in the numerical 
implementation of Eqs. ( 1)— (6). This is because the 
integrands are rapidly varying functions of 0 and 4>> 
i.e., because of the P,i m i(cos 0) and exp(-im<|>) factors 
and the fact that and are proportional to 

exp(i£r cos 0) with the usual evaluation criterion r — 
a. Thus dense grids are required in both the 0 and 
the 4) directions to obtain convergence of the numeri- 
cal integrations. 18 The localized model approxi- 
mates the integrals in Eqs. (5) with an analytical 
function, thus decreasing the computer run time 
many fold. For an off- axis- focused Gaussian beam, 
the localized model of the beam-shape coefficients is 17 


The accuracy of Eqs. (7)— (12) in approximating Eqs. 
(5) for a focused off-axis Gaussian beam was exam- 
ined in Refs. 15 and 17-19. Equations (7)— (12), 
however, are not in optimal form for numerical 
computations. In particular the Kronecker delta 
implies that for a given value of j, only one value ofp 
is considered. When this is taken into account, the 
resulting infinite series in j is recognized as that of a 
Bessel function. 20 The localized beam-shape coeffi- 
cients A, m loc and B [m [oc may then be written in the 
more compact form 



r A, +loc 
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Similarly, 


B b loc 


*lm 


(B, m +loc if m > 0 
B tm - ]oc ifm < 0 
B ;o loc if m ■ 0, 


(17) 


where 


B rloc — 
Im “ 


±Fi j 




(«/ m _i(P) + (r 


s,. 1 ” - F 


i.l) v*’?* 


(18) 


Equations (14H18) simplify further when the beam 
waist and the particle orientation are such that x,- ^ 0, 
Vf = 0 or when Xf — 0, * 0 [see Figs. 1(b) and 1(c)]. 

If Xf = Vf = 0, the coefficients reduce to their values in 
the on-axis localized model. 16 ^ 
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Because Bessel functions and modified Bessel func- 
tions are related to each other by 21 

J n (x) = i n I n (-ix\ (19) 

Eqs. (14H18) may be rewritten as 

\ m— 1 

[I m -i(Q) + (r, s ) 2 7 m+1 (Q)], 
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**-lm 1 l 


-irp 
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where 


Q = 2s\ l + g 


l\(x f 2 + y f 2 Y /2 l 2isz f y i 


Wo 


1 - 


w 0 


= -iP. (22) 


Again Eqs. (20) and (21) further simplify when x f * 0, 
y f = 0 or when x f = 0 ,y>f* 0. Because these special 
cases are examined at length in Sections 3 and 5, 
below, we present the simplified expressions here. 
If Xf * 0 and y f = 0, Eqs. (20) and (21) become 


A, m £loc = Vi 


— I x 
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B/o'“ - 0, 

and ifx/ = 0 and y f * 0, Eqs. (20) and (21) become 
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In Section 3 we find Eqs. (14), (16), and (18) with 
Bessel functions of a complex argument to be the 
most efficient form of the localized model coefficients 
when the Gaussian beam has spread beyond its focal 
waist at the position of the particle. Similarly we 
find Eqs. (20)— (24) with modified Bessel functions of a 
complex argument to be the most efficient form when 
the particle is within the focal waist of the beam. 

3. Numerical Implementation of the Localized Model 
for Gaussian Beam Scattering 

In Section 3 we discuss three numerical aspects of the 
computation of the far-field scattered intensity of 
Eqs. (1) and (2). These are (a) the evaluation o{J n (P) 
and 7„(Q) in the expressions for the beam-shape 
coefficients; (b) an algorithm patterned after Wis- 
combe’s method in plane-wave Mie theory 3 for calcu- 
lation of the angular functions ir/ |m| (0) and T/i m i(0); 
and (c) the finding that the entire range of m values in 
Eqs. (2), i.e., -I <> m <, l, need not be computed. 
We determine the value of for which inclusion of 

the terms —m mia < m < m max in Eqs. (2) yields 
accuracy of 1 in 10 8 in the computation of 7(r, 8, 6). 


A. Modified Bessel Functions 

Consider the modified Bessel function 7„(Q), where n 
is a nonnegative integer and Q is complex. Refer- 
ence 22 states that for Re Q < 12 or Re Q < n, the 
Taylor series expansion 


UQ) = 


(QY ^ (Q 2 /4)* 

\2j j~ 0 k\(n + k)l 


(25) 


is rapidly convergent, and for Re Q > 12 and Re Q > 
n, the asymptotic series 


7„(Q) = 


exp(Q) 

(21I-Q) 1 ' 2 


1 + 2 


(- 1 )* 

A!(8 Q) k 


(4n 2 — 1) 


x (4 n 2 - 9) . . . [4n 2 - (2k - l) 2 ] 


(26) 


may be used efficiently. No upper limit for the k sum 
is given in Eq. (26) because the asymptotic series 
diverges; i.e., the terms of the series become smaller 
and smaller, reach a minimum size, and then become 
larger and larger. Using the remainder theorem for 
an alternating series, 23 we achieve the best approxima- 
tion to I n (Q) when all the terms up to and including 
one before the smallest term are summed. Examin- 
ing Eqs. (25) and (26), we found that 1 in 10 8 
convergence could be achieved only if Arfken’s re- 
gions of applicability were changed to Re Q < 12 or 
Re Q <; n + 2 for Eq. (25) and Re Q > 12 and Re Q > 
n + 2 for Eq. (26). The upper limit of the k sum in 
Eq. (25) for 1 in 10 8 convergence was found to be 

k max = Re Q + 7 + 0.5 1 Im Q | . (27) 


M 


For Eq. (26) the upper limit of the k sum for the same 
convergence criterion was found to be 

- » + 12. (28) 


B. Bessel Functions 

Consider the Bessel function J n (P). If Re P ^ 12 or 
Re P ^ n + 2, the Taylor series expansion 


For Im Q = 0 we checked our computed values of 
I n (Q) by comparing them with the tabulated values in 
Ref. 24. For Im Q * 0, the computed values of I n {Q) 
were checked to make sure that they satisfied the 
symmetry relation 25 

/„(<?*) = (29) 

The value of in Eqs. (27) and (28) is valid only if 
Re Q > \ Im Q | , or when 



{ — 1 ) a (jP 2 / 4)* 
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is rapidly convergent. Convergence of 1 in 10 8 was 
achieved when the upper limit of the k sum was 

4 max = ReP + 12 + 0.5 ImP (32) 
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where L is the spreading length of the Gaussian 
beam. For Re Q < | Im Q | it was found that the 
number of terms required for convergence rapidly 
increased, rendering this method of computation 
inefficient. 

Within the beam focal waist given by Eq. (30), if 
z f > 0, the beam is converging at the position of the 
spherical particle and Im Q > 0. If z,- = 0 , the beam 
focuses in the plane that contains the particle, and 
Im Q = 0. If Zf < 0, the beam is diverging at the 
position of the particle, and Im Q < 0. The position- 
ing of the particle in the beam and its corresponding 
location in the complex Q plane are illustrated in 
Fig. 2. 


or | Re P| > Im P. For Re P > 12 and Re P > n + 2, 
the asymptotic series is 26 


Jn(P) = 
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Fig. 2. (a) Complex Q plane as defined in Eq. (22). (b) Focused 

Gaussian beam that is propagating from left to right. The points 
labeled A through E in the off-axis beam are the positions of the 
spherical particle within the beam's focal waist. They also corre- 
spond to the indicated locations in the complex Q plane. 


(4 n 2 - l)(4n 2 - 9)(4n 2 - 25) 


3!(8P) 3 


x sin 

H)* 


L 2 J 



(34) 


We call the first term in each series in Eq. (34) the k — 
— 1 term, the second term k = 1, the third term k = 3, 
etc. Equation (34) also diverges if too many terms 
are considered. But when |ReP| > ImP, it con- 
verged to 1 in 10 8 for 


*max 33 n + 9 * 

Again, when | Re P| < Im P, the number of terms in 
Eqs. (31) and (34) required for 1 in 10 8 convergence 
grew rapidly, rendering this method of computation 
inefficient. 

For realP, we checked our computed values of J n (P) 
by comparing them with the tabulated values in Ref 
27. If P is complex, care must be taken in the 
evaluation of Eq. (34). For the trigonometric func- 
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tions in Eq. (34), we have 


ImP 


COS 


P - 




= COS 


ReP- 



cosh(Im P) 


- i sin 


ReP- 



sinh(Im P), 


sin 


P - 



2 


= sin 


Re P - 



cosh(ImP) 


+ i cos 


ReP - 



sinh(Im P). 


(36) 


But the P 1/2 factor in the denominator is potentially 
problematical because the square root of a complex 
variable is a two- sheet function. Thus the question 
arises as to which sheet a particular value of P is on. 
If the beam focuses upstream from the particle with 
Zf < 0, Eq. (16) yields Re P > ImP and Im P > 0. 
In this region of the complex P plane adjacent to the 
positive real axis we have 

P = re* (37) 


(a) 



Fig. 3. (a) Complex P plane as defined in Eq. (16). fb) Focused 

Gaussian beam that is propagating from left to right. The points 
labeled A through D in the off-axis beam are the positions of the 
spherical particle outside the beam’s focal waist. A-D also corre- 
spond to the indicated locations in the complex P plane. 


C. Angular Functions 

The angular functions T7,i m i(0) and T,' m| (9) satisfy the 
recursion relations 18 


21 + 1 


*i + r<A) = z + 1 


cos 0'T7/ m (0) 


l + m 


6 ), 


(41) 


l + 1 - m 

T, m (0) = i cos e-, m W - (t + TO)-,_r(e), (42) 

with the starting values 


pi /2 = r i /2 exp(i0/2). (38) 

On the other hand, if the beam focuses downstream 
from the particle with Zf > 0, Eq. (16) yields —Re P > 
Im P and ImP > 0. The positioning of the particle 
in the beam and its corresponding location in the 
complex P plane are illustrated in Fig. 3. Contrary 
to the situation in Fig. 2, the z f > 0 region is disjoint 
from the Zf < 0 region, but satisfies the symmetry 
relation 

(P) Zf >o = ~(P%<o- ( 39 ) 

Thus, rather than computing J n {P) separately for 
Zf > 0 in the disjoint region, we compute it by the use 
of the identity 28 

Jn[{P) tf >o} = (-U n Jn*[(P), f <ol (40) 


-r7,_i'(0) = 0, 

it/( 0) = (2/ - 1)!! sin /_1 0 (43) 

for m > 1. For m = 0, the starting values are given 
in Ref. 18. Often one wishes to calculate the far-field 
intensity at many angles 0 to construct an angular- 
scattering diagram. Because the computation of 
Tr; m (0) and T/ m (0) is within a triple do loop, i.e., m y l , 
and 0, savings in computer run time may be realized if 
the number of multiplications within the triple do 
loop is minimized. Following Wiscombe, 3 for m > 1 
we compute the angular functions with 


s = cos 0n / m (0), 

(44) 

t = s- n,_r(e), 

(45) 

i 

t,™( 0) = -t- n,_r(ei, 

m 

(46) 
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n i+ r(0) = s + 


l + m 
l + 1 - m 


T, 


(47) 


where 


n, m (0) = ( 48 > 

If the values of l/m and l + m/(l + 1 — m) are 
precalculated, we can use Eqs. (44)— (47) to compute 
n /+l l'»i( 0 ) and T^ m, ( 0 ) with 3 multiplications and 3 
additions, whereas Eqs. (41) and (42) require 6 multi- 
plications and 2 additions. No decrease in the num- 
ber of multiplications is obtained, however, by compu- 
tation of the quantities Si + S 2 and S x - S 2 rather 
than computation of Si and So individually as in Eqs. 
(49), below. 

D. Number of Terms in the m Sum 

In previous Gaussian beam-scattering calculations, 18 * 29 
it has been useful to interchange the order of the l and 
m sums in Eqs. (2), yielding 


*<•.♦> -fsF7 ) fWWW 

mmax /max 0/ 4* 1 

+ 2 2|rTTi-n,-(»» 

X [-A,„ + exp(mi<f>) + A tm - exp(-tm4))j 

^mu /max 0/4- 1 

x [5 Im + exp(im<i>) + B im ~ exp(-:md>)] 


of l, A lm ~ and B lm ~ rapidly decrease, and 11, '"(9) and 
T ( m (0) rapidly increase, but the product of the beam- 
shape coefficients and the angular functions, which 
we call the weighted beam-shape coefficients, also 
rapidly decreases. 18 This permits truncation of the 
m sums in Eqs. (49) at m maj( « l with little loss in 
accuracy. 

This result may be demonstrated analytically as 
follows. Consider for simplicity the case z f = 0 so 
that Q is real. Consider also a relatively high partial 
wave, a small beam focal waist, or a relatively large 
off-center impact of the beam on the sphere so that 
Q > 12. If x f * 0 and y f = 0, the beam-shape 
coefficients are given by Eqs. (23). We now examine 
these expressions as a function of m. When m is 
small, / m *i(Q) may be approximated by the first two 
terms of Eq. (26), yielding 


|A im = l0C | 
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x (2u-Q)" l/2 exp(Q)^2 - Jq ~ ~qJ ' 


exp 

[-44)1 

exp( -x f 2 / w 0 2 ) 

/ i' 

\ l + 2 

\ A 

jm- l 


x (2irQ) _1/2 exp(Q)|-Q-j • (52) 


/ max 21 + 1 

S 2 (9> 4>) = 2 ^TT) A,oa,T,0(9) 

mmj, /m« 21+1 




x [A tm * exp(tm<b) + A /m exp(-i/mf>)] 

^max /max 9/4-1 

x [Bim* exp(im<i>) - B lm ~ exp(-im<J>)], 

(49) 

where the largest partial wave / max is given by 3 - 4 

/ max = 2 + X + 4.3X 1/3 , (50) 

the size parameter X is 

2i7a 


X = 


(51) 


and the value of m max is yet to be determined. For 
the examination of low-order MDR s, the value of l max 
may have to be increased somewhat. 3 Numerically 
it has been found that as m increases for a fixed value 


Similar expressions occur if x f = 0 and 0. For 9 
away from 0° or 180° and for large /, the angular 
functions approach the asymptotic values 18 

n, m (0) - / m+l/2 (sin 9)" 3/2 



7 1\ m— tt 

X COS 

[i + 2) 0 + T ~ 4 


/ 2 \l/2 

T/ m (0) (-1 Z m+l/2 (sin 0)' l/2 



7 1\ VI TT TT 

x sin 

\ + 2/ + ~2 4 


or, re tain ing only the powers of l in expressions ( 53), 

l7/*(e) i = ^inrwi = ( 54 ) 

We can now see from expressions (52) that the 
beam-shape coefficients decrease as a function of m as 
whereas from expression (54) the angular 
functions increase as Z m+l/2 . Thus for small m the 
dominant m dependence of the two contributions 
cancels, and the weighted beam-shape coefficient 
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\A lm ^-r I of Eqs. (49) decreases only quadratically in 
m. Likewise |B (m - loc T/ m | starts out much smaller and 
increases linearly in m. Because from expression 
(54), rriT:i m <k x, m for small m, the weighted beam- 
shape coefficients |A; m =loc m-n- ( m | and |B; m rloc m-rr ( m | are 
much smaller than |A; m - loc -r ( m | and |B; m =loc T; m | . 

When m is large, may be approximated with 

Eq. (25) as 

exp -s 2 ll + -j exp( -Xf 2 /w a 2 ) 
|A, m =‘«| - |B (m " l0C | = 


IQ /9} m ~ l 

X- — exp(Q 2 /4m). (55) 

(m - 1)! 



shape coefficients are roughly equal and decrease 
almost exponentially as a function of m. 

Figure 4 suggests that the m sum in Eqs. (49) may 
be cut off for large Q and l at m raax <sc l with minimal 
loss in accuracy. The cutoff value rn max was deter- 
mined by the criterion 


Mm ~l lamallm 


J— 7TJ exp(Q 2 /4m max ) 

(m max - 1)! 

x (2-irQ) 1/2 exp(-Q) < 10 -8 . (56) 


The solution of Eq. (56) was obtained numerically, 
and it is well approximated by the relation 


m max = 6.5 Q 1 ' 2 for 6 < Q < 40. (57) 


A similar expression occurs if Xf = 0 and yy * 0. Thus 
the weighted beam-shape coefficients are approxi- 
mately equal and rapidly decrease when Q 2 a 4m. 
To check this predicted dependence of the weighted 
beam-shape coefficients on m, we computed |A (m - loc r/ m | 
and |J5 /m = loe T, m | numerically with Eqs. (23), (25H28), 
and (48) for l = 430, x f = 42.5 n.m, y f = z f = 0, w 0 = 
13.3 n-m, and X = 0.6328 M-in corresponding to Q - 
20.83. The results are shown in Fig. 4. As pre- 
dicted in expressions (52), initially |A im rloc T( m | is large 
and slowly decreasing, and |B/ m - loc T( m | is small and 
increasing. For m > 14 the two weighted beam- 



Fig. 4. Weighted beam-shape coefficients |A, m i “ : T,' n | (filled circles) 
and | B /m \ ( open circles i as a function of | m j for / = 430 and an 
off-axis Gaussian beam with k = 0.6328. 


This result was tested with Eqs. (20)— (22) and (25)— 
(28) and expression (54) to compute the weighted 
beam-shape coefficients. Equation (57) was found to 
be accurate in every instance. When Q was made 
complex, the weighted beam-shape coefficients fell by 
a factor of 10® when 

rn max = (6.5 + 2.0^|^)(ReQ) 1 '' 2 

for|ImQ| ^ Re Q. (58) 

As a final check of Eqs. (25H28), (31) and (32), 
(34)-(40), and (58) we calculated the far-field scat- 
tered intensity for z f = w 0 /2s, the boundary between 
the modified Bessel function representation and the 
Bessel function representation of the beam-shape 
coefficients. The far-field scattered intensity was 
calculated with each representation, and the results 
agreed to better than 1 in 10 8 . 

E. Inclusion of the Incident Beam 

Before we examine the run time of our Gaussian 
beam-scattering computer program, an important 
addition to the scattering amplitudes of Eqs. (2) and 
(49) must be made. Equations (2) and (49) give us 
the amplitude for the far-field scattered portion of the 
electromagnetic fields exterior to the spherical particle. 
But the entire exterior field, the scattered field plus 
the incident field, is measured in the experiments. 
Thus the amplitude of the incident field should be 
appended to Eqs. (2) to have the resultant expression 
agree with experimental observations. The incident 
beam is included in far-field plane-wave Mie theory 
only at 0 = 0 but is in included in the near-forward 
direction in the near field. 30 For Gaussian beam 
scattering the incident beam must be included in both 
the near field and the far field in the near-forward 
direction because of the spreading of the incident 
beam. 3132 

In particular consider a narrow Gaussian beam 
that is incident only slightly off-axis upon a large 
particle. Because the particle obstructs most of the 
incident beam, the near-forward diffracted field should 
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be quite weak because only the tail of the Gaussian 
beam passes the edge of the particle. Yet the 
numerical implementation of Eqs. (1). (2), and ( ) 

yields a large unobserved diffraction like peak m the 
near-forward direction, 32 which is canceled by the 
spreading of the incident beam. The addition of the 
incident beam to the scattering amplitudes for an 
off-axis beam is given by 

for 0 > 10 s 


S 1 totai (e, <b) = 


S 2 total (0, *) = 


5 ^ scattered/ Q ; 

s 1 scat t«ed(e i 4>) 

- sin <i> S inadent (9. <t>) for 0 ^ 10 s, 

[ S 2 scatten,d (0, <j>) f° r 9 > 10 s 

S 2 ’ cattered (0, 4>) 

- cos d> S in<:idem (0, <i>) for 0 ^ 10 s, 


ginddent(0 > <|>) = — exp( 


x exp 


-0 2 /4s 2 ) 


i0 


x f yf . , 

— cos <b + — sin cp 

,W 0 W 0 


x exp(— iZf/swo)exp(id 2 Zf/2sw 0 ). (59) 

We obtain the corresponding equations for an on-axis 
beam from Eqs. (59) by setting x f = y f - 0. The 
inclusion of the incident beam is implemented m our 
Gaussian beam- scattering computer program. 

4. Timing Study of the Localized Model of Gaussian 
Beam Scattering 

To assess the performance of our Gaussian beam- 
scattering program, we tested it on the situation in 
which a = 50 pun, n = 1-333, X = 0.6328 p.m w 0 - 10 
p.m, x f = z f = 0 ,y f = 20 pm, 4> = 90“ and for 361 values 
of 0 in the interval -180“ ^05 180 . The size 
parameter for this case is X - 496.46, the largest 
partial wave is = 532, the degree of beam 
confinement is s = 0.01, and the off-centeredness of 

the beam is = 21.3. According to Eq. (57) this 
value of corresponds to m mex = 30, vduch was 
used as the upper limit of the m sum m Eqs. (49). 
This program, as well as all the other programs for 
which timing studies were made, was run on a 
Compaq 386-33 MHz personal computer equipped 
with a Weitek numerical coprocessor. The run time 
of the localized model Gaussian beam program was 
195 s for the parameters given above. Less than 1 s 
of this time was spent computing the incident beam of 
Eqs (59) and the Mie partial-wave scattering ampli- 
tudes for 1 ^ l Z /max- Because the localized approxi- 
mation replaces the numerical integrations or Eqs. 
(51 only 9 s were spent computing A| m =loc and B tm - w 
for 1 < l ^ Imox and 0 £ rn <, with Eqs. (24). 
The program spent 3:05, or almost 95% of the run 
time, computing irtfAJ and ^(0), multiplying the 
beam-shape coefficients by the angular functions, and 
adding everything together to obtain the scattering 
amplitudes. This division of run time is similar to 


that reported by Wiscombe for plane-wave Mie theory. 3 
When z f * 0, the beam-shape coefficients were calcu- 
lated with Eqs. (20) and (21) rather than the simpler 
Eqs. (24). The fact that Q was complex meant that 
20 s more was required for the calculation of A( m 
and B, m - loc for 1 £ l <. l max and 0 < m < m max . The 
time for all the other computations was unchanged. 
For comparative purposes, a plane-wave Mie theory 
calculation for a = 50 |im, n = 1.333, X = 0.6328 M-m 
and for 361 values of 0 in the interval -180“ ^ 9 ^ 
180“ took slightly less than 3 s on the same computer. 
Thus our Gaussian beam program runs almost 70 
times slower than Mie theory for these parameters. 

In Eqs. (2), the full range of l and m values is 1 £ 
l ^ and -l < m <, l. For the parameters of our 
numerical experiment this requires the computation 
and storage of 568,178 beam-shape coefficients. 
Truncating the m sum at m max = 30 reduces the 
number to 63,166, which is 11.1% of the total. 
Using the symmetry-relations for Ai m loc and Ai m 
and for and B (m ‘ loc of Eqs. (14), (18), (20), and 

(21) further reduces the number of coefficients com- 
puted to 32,116, which is 5.6% of the total. Thus the 
truncation of the m sum at 1 in 10 8 accuracy for the 
far-field intensity represents a substantial savings in 
computer run time in this example. 

A program that also computes the scattering ampli- 
tudes with Eqs. (49) but computes Ai m and B lm with 
numerical integration of Eqs. (5) was written. The 
grid size for the 0 and <j> integrations required for 
convergence of the numerical integrations is given 
elsewhere. 18 For the parameters of our numerical 
experiment with = 30, the run time for this 

program when 32,116 beam-shape coefficients were 
used was 4.5 h which is a factor of 83 slower than our 
localized model program. If the full range of m 
values had been used, the run time would have been 
longer by another factor of at least 17.8. 

In Refs. 9-11 the incident Gaussian beam is ex- 
panded in an angular spectrum of plane waves. The 
plane waves are then decomposed into vector spheri- 
cal harmonics. We obtain the total vector spherical 
harmonic coefficients by summing the individual 
plane-wave coefficients over the angular spectrum. 
The total vector spherical harmonic coefficients are 
then input into a T-matrix program for calculation of 
either the far-field intensity or the interior source 
function. 29 The computer run time required for the 
computation of each set of total vector spherical 
harmonic coefficients a €m n, a omn> and b omn was 

1.57 s. Thus computation of the 32,116 sets ot 
coefficients required for our test situation takes 14 h. 

The Rouen computer program for Gaussian beam 
scattering described in Ref. 19 is in many aspects 
s imilar to the program that we have described here. 
The localized approximation is used and Eqs. (7) are 
written as a single sum over j. But (a) the sum was 
not recognized as being a Bessel function or a modi- 
fied Bessel function, (b) the series was truncated 
when an individual term fell below 10' 30 rather than 
at 1 in 10 8 accuracy, and (c) the number of m values 
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was set at m max = 10 rather than at m max of Eqs. 
(57}-{58). This resulted in A lm ]oc and B im loc being 
computed to much greater than 1 in 10 8 accuracy. 
But for Q > 2.4 a number of weighted beam-shape 
coefficients were omitted that were larger than 10 “ 8 
of the coefficients that were included. For the param- 
eters of our numerical experiment, the Rouen pro- 
gram computed the localized approximation beam- 
shape coefficients in 96 s which is a factor of 10.7 
slower than with our program. But, when the local- 
ization approximation is used, the most time-consum- 
ing part of the program is the computation of T/ m , and 
T7/ m . Thus the entire Rouen program is only a factor 
of 1.45 slower than ours. 

The run-time study described here was for only one 
particular example of Gaussian beam scattering, and 
it would be unwarranted to extrapolate the compari- 
son between the various computational schemes to all 
cases of Gaussian beam scattering without further 
testing. F or example, consider the published calcula- 
tions of Gaussian beam scattering given in Table 1. 
Of particular interest in Table 1 are the values of Z^, 
the highest partial wave in the computation, and 
Qmax> which is a measure of the highest partial wave, 
the extent of beam focusing, and the degree of 
ofF-centeredness of the incident beam. In each of the 
references cited in Table 1, the authors had a differ- 
ent goal in mind when performing the calculations, 
thus dictating different choices for Z max and Q max . 
Lock 18 was interested in rainbow formation and thus 
required a large particle, i.e., Z max = 565. Because 
the particle was large, the incident beam did not need 
to be tightly focused, i.e., Q ^ = 9, to see the effect for 
which he was looking. The combination of large Z max 
and relatively small is tailor-made for both the 
localized model and the truncation of the m sum 
because computing the full range of A im and B lm 
coefficients by numerical integration for large l would 
take a prohibitively long time on a personal computer. 
The same consideration holds true in computation of 
calibration curves for particle-sizing instruments in 
the large-particle regime. 36 

Barton et aZ. 8 * 33 * 34 were interested in smaller par- 
ticles, i.e., Z max = 45. But for them to see the effects 
that they were looking for, the incident beam had to 
be more tightly focused, i.e., <J max = 37. Their beam 
localization of s = 0.084 is near the limit of the 
validity of the localized model for an off-axis beam, so 


Table 1 . Parameters in Published Gaussian 
Beam-Scattering Calculations 


Reference 

18 

a (p.m) 
43.3 

\ (p-m ) 
0.5145 

u>o(^m) 

20.0 

s 

0.004 

^max 

565 

Qmax 

9.0 

8 

2.5 

1.06 

2.0 

0.084 

27 

5.8 

33 

5.0 

1.06 

2.0 

0.084 

44 

37.4 

34 

2.5 

0.5145 

1.0 

0.082 

45 

37.3 

10 

8.0 

1.06 

2.0 

0.084 

64 

86.5 

11 

8.0 

1.06 

2.0 

0.084 

64 

64.9 


13.4 




100 

113.4 

35 

9.45 

1.06 

2.0 

0.084 

74 

88.7 

19 

38.08 

0.5145 

10.0 

0.0082 

500 

8.2 


its use still yields a great decrease in computer run 
time. But now because m max = 40, there is not 
much point in truncating the m sum before m = 1. 
Barton et al. performed their computations on a 
Silicon Graphics 4D/380 VGX Super Computer, 37 
taking advantage of the much faster speed of 34 
MFlops to perform the numerical integrations for A lm 
and B[ m in Eqs. (5). Khaled et aZ. 10 - 11 - 35 were also 
interested in small particles, i.e., Z^ = 75, and in 
particular with MDR’s excited by a tightly focused 
beam, s = 0.084, that was incident upon a spherical 
particle somewhat beyond the particle’s edge with 
Qmax ~ 100. Again the s value is just within the 
range of applicability of the localized model, but 
truncation of the m sum again is not possible. If the 
beam had been even more tightly focused, with s > 
0.1, the localized beam model might produce a rela- 
tively poor approximation to the actual beam pro- 
file. 15 

The moral of the story is that the computational 
cost of Gaussian beam-scattering calculations de- 
pends on the size of the spherical particle through 
Z max , on the degree of focusing and the degree of 
off-centeredness of the incident beam through Q^, 
and on the speed of one’s computer. The computa- 
tional method described here provides the greatest 
computational savings for Q max < 40 and large Z^, 
thereby permitting the computation to be easily 
handled by a personal computer. For much larger 
values of Q ^ or smaller values of l maX) the computa- 
tional savings may not be as great. But the calcula- 
tion is still very efficient on a personal computer. 

5. MDR Excitation in the Localized Model 

It has been found both experimentally 38 * 40 and theo- 
retically 10 * 33 within the last few years that one can 
greatly enhance the MDR excitation rate in a dielec- 
tric microparticle by having a tightly focused Gauss- 
ian beam that is incident somewhat beyond the 
microparticle’s edge. Theoretically the reason for 
this is that the spherical Bessel function ji{nkr) inside 
the particle couples to the spherical Bessel function 
ji(kr) outside the particle. 10 Because the energy 
density of a MDR is largest just inside the particle 
surface, we must have Z > nka to make ji(nkr) reach 
its peak value there. 41 But because the classical 
impact parameter b of the geometric light ray associ- 
ated with the particle wave Z is 42 * 43 l ~ kb, we needed 
b/a < n for the excitation of the MDR. This impact 
parameter describes a light ray that passes the par- 
ticle somewhat beyond its edge. 

The localized model permits a more detailed descrip- 
tion of the effect. Consider a focused Gaussian beam 
with Xf * 0 and yf - Zf = 0. Assume that the 
particle’s radius and refractive index and the beam’s 
wavelength are such that the Mie particle-wave scat- 
tering amplitude a/ resonates, producing a TM-type 
MDR. In Eqs. (4) the effect that a t has on the 
far-field intensity is modulated by the beam-shape 
coefficients A lm for -Z < m < L If the A lm are small, 
the MDR is suppressed. If the A lm are large, the 
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strength of the MDR is enhanced. The same holds 
true for the interior source function because the Mie 
interior amplitude c* and the scattering amplitude a< 
resonate simultaneously and because the interior 
field is proportional to ciAi m . As can be seen in Eqs. 
(23) and Fig. 4, when x f * 0 and yy = 0, we have 
^ m loc » Bf m loc for small m. Thus the TM resonances 
are orders of magnitude stronger than the TE reso- 
nances. On the other hand, when x, = 0 and y f * 0, 
Eqs. (24) show that B [m ]oc » A lm ioc for small m. 
Thus the TE resonances that are proportional to 
B lm loc bi are orders of magnitude stronger than the TM 
resonances. 

Consider the weighted beam-shape coefficient 
=loc T; m of Eqs. (23) and expression (54). When m 
is small and m.ay be approximated by the first 

two terms of Eq. (26), we have 


1 

x f 

\-l/2 

f 

[ / i) 

*r 

12' 

=5 /( 4tts 

Wq 

j exp' 

l" 

H' + 5) - 

w 0 




3 m 2 \ 

4Q" O']' 


(60) 


As seen in Fig. 4, Eq. (60) decreases quadratically as a 
function of m. But as a function of x f Eq. (60) peaks at 



|x f | 2 t r 


(61) 


Because the resonating partial wave l for a low MDR 
order number i occurs for 44 - 46 


** = (z + l) + o,(/ + ^) 1/3 2- 1/3 - V(n 2 - l)- l/2 , 
n for a TE resonance 

V= 1 , _ . , (62) 

- for a TM resonance 

n 

where 46 Ai(-ai) = 0, the greatest enhancement in the 
strength of the MDR occurs for an incident beam with 

y f fay / 3 1 n 1 /2\ 1/3 

a = n ~ ai \2) X* /3 + (n 2 - l) l/2 X + 6\n) F' 2 

for a TE resonance 


Xf ln\v* 1 1 l/2W3«.-» 

a " n ~ ai \2j X 2 ' 3 + n(n 2 - 1) 1/2 X 6\n) X 

for a TM resonance, (63) 

or for focusing somewhat outside the edge of the particle. 
As an example, the TEsa,! resonance has 
L = 58, i = 1, and X - 47.3094299 for a = 1.36. 
Expressions (63] yield y f /a = 1.24, which agrees wefl with 
the value of 1.23 in figure 2b of Ref. 10. Similarly the 
TM341 resonance has 7 = 34, t = 1, andX = 29.753 for 
m = 1.33. Expressions (63) yield x f /a = 1.16, which 
agrees well with the value of approximately 1. 14 in figure 
7 of Ref. 33. Also the TE^i resonance has l = 34, i = 1, 


and X = 29.365 for m = 1.33. Expressions (63) yield 
y f /a = 1.18, which agrees well with the value of approxi- 
mately 1.18 in figure 8 of Ref. 33. Expression (60) also 
shows that the low azimuthal modes m of the MDR are 
all excited to nearly the same extent, whereas the degree 
of excitation of the high azimuthal modes fall off rapidly. 
This is illustrated in Fig. 5. 

In summary, the most important result of this 
paper is that Gouesbet’s localized model for the 
beam-shape coefficients in scattering of a focused 
Gaussian beam by a spherical particle may be written 
in terms of either Bessel functions or modified Bessel 
functions. On the one hand this simplified form 
leads to the construction of a fast-running computer 
program for Gaussian beam scattering that can be 
implemented on a personal computer. On the other 
hand the simplified form provides a simple analytical 
formula for the beam-shape coefficients that permits 
one to obtain an intuition of various effects that occur 
in Gaussian beam scattering. 
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Fig. 5. Weighted beam-shape coefficient as a function 

of Xf for l * 430 and an off-axis Gaussian beam with X = 0.6328 M.m 
and y f = Zf = 0. The short dashed curve is the weighted 
beam-shape coefficient for jm| = 1. This is the only coefficient 
that remains nonzero in the on-axis limit and is proportional to the 
MDR excitation rate by an incident plane wave. 
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